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We study the stability of unstable steady states in scalar retarded time-delayed systems subjected 
to a variable-delay feedback control. The important aspect of such a control problem is that time- 
delayed systems are already infinite-dimensional before the delayed feedback control is turned on. 
When the frequency of the modulation is large compared to the system's dynamics, the analytic 
approach consists of relating the stability properties of the resulting variable-delay system with 
those of an analogous distributed delay system. Otherwise, the stability domains are obtained by a 
numerical integration of the linearized variable-delay system. The analysis shows that the control 
domains are significantly larger than those in the usual time-delayed feedback control, and that the 
complexity of the domain structure depends on the form and the frequency of the delay modulation. 
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I. INTRODUCTION 

In spite of the fact that the control problems have been 
thoroughly investigated from a theoretical aspect and the 
results are being iniplemented in concrete real systems 
for several decades [l| , the control of chaotic dynamical 
systems is a relatively new area of research. The appear- 
ance of the pioneering paper by Ott, Grebogy and Yorke 
(OGY) in 1990 boosted quite an interest among nonlinear 
scientists [2|, being a reason for a large number of pub- 
lished papers on chaos control f^-Q- The OGY method 
utilizes the existence of infinitely many unstable periodic 
orbits (UPO) within the structure of the chaotic attrac- 
tor, applying a small externally controlled perturbation 
to suitably chosen parameters of the system when the 
trajectory is in the neighborhood of an UPO whose con- 
trol is desirable. The system is then externally forced to 
follow otherwise unstable behavior corresponding to that 
UPO. The numerical simulations and the experimental 
implementations showed that the method by itself has 
some drawbacks concerning the robustness with respect 
to the external noise and its practical realization, since 
it requires a continuous monitoring of the evolution of 
the system from the outside and the knowledge of the 
equations that describe the system's dynamics. 

The OGY idea stimulated a development of a rich va- 
riety of new chaos control techniques. Among those is 
the time-delayed feedback control (TDFC) proposed by 
Pyragas in 1992 [3,l3i shown to be much more fiexible for 
practical purposes with respect to OGY (the monitoring 
of the system and the knowledge of the exact positions 
of UPOs are not required) and quite robust against the 
effects of noise. The control force is applied as a con- 
tinuous feedback proportional to the difference between 
the current state of the system and the state of the sys- 
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tem delayed by the constant time T. If the time delay T 
is chosen to coincide with an integer multiple of the pe- 
riod of the target UPO, then the control force will vanish 
when the target state is reached and the control is nonin- 
vasive. For stabilization of unstable steady states (USS), 
the choice of the parameter T is not as restrictive as in 
the case of UPOs, and the interval of T for which TDFC 
is successful is shown to be system-dependent |9l-[ll|. 

In parallel to various practical applications of TDFC 
[l2 - [2l| . an effort has been put into progress to gener- 
alize or modify the original control scheme in order to 
improve its performance. Some extended TDFC schemes 
employ multiple time-delays to stabilize strongly unsta- 
ble periodic orbits |22 h25[ . Others are introducing an 
oscillating feedback gain [26[ or an extra unstable degree 
of freedom in the feedback loop 271-2911 to overcome the 
so-called odd- number limitation 30l-l32|j which was re- 
futed recently [33l - l36| . In a recent work [33], it has been 
shown that the efficiency of the TDFC method to con- 
trol USS can be significantly improved by including a 
variable time-delay into the TDFC scheme in a form of a 
deterministic or stochastic modulation in a fixed interval 
around a nominal delay value. Stochastic changes in the 
delay time are natural due to the omnipresent noise in 
any physical system. In the circumstances, the enhance- 
ment of noise along the delay line could be desirable as 
it is leading to improved stability of the system. On the 
other hand the modulated delay described by some peri- 
odic function could be realized by periodically changing 
some characteristic distances in electric or laser systems 
by introducing piezoelements. This variable delay feed- 
back control (VDFC) has been shown successful in stabi- 
lization of USS in low-dimensional chaotic systems using 
different types of delay modulations. The ongoing analy- 
sis shows that VDFC can also improve the control domain 
of UPOs with respect to TDFC for a specific choice of 
the delay modulation j38l |. 

The purpose of this paper is to investigate the effects 
of stabilization of unstable equilibria by a variable-delay 



feedback control in a class of nonlinear dynamical sys- 
tems described by scalar retarded delay-differential equa- 
tions (RDDE) involving the value of the state variable 
at a discrete time lag. A delay differential equation is 
called retarded if the highest order derivative only occurs 
with one value of the argument, and this argument is not 
less than the arguments of the unknown function and its 
lower order derivatives appearing in the equation 39 41]. 
In contrast to low dimensional dynamical systems, delay 
differential equations are infinite dimensional, since it is 
necessary to specify a continuum of initial conditions over 
the interval length equal to the time delay. The interest 
for such equations is caused by their frequent occurrence 
in numerous physical, biological and engineering models, 
where the time delays are a natural manifestation of the 
system's dynamics [42h45| | . 

The paper is organized as follows. In Section II we 
perforin a linear stability analysis of USS in the free run- 
ning RDDE system and in the system under VDFC. The 
frequency of the delay modulation in the feedback loop 
is considered to be sufficiently large compared to the in- 
trinsic timescale of the unperturbed system, allowing an 
approximation of the variable delay system with a dis- 
tributed delay system 46|. In Section III, we numerically 
illustrate the VDFC method in the chaotic Mackey-Glass 
system. The domains of successful control are first com- 
puted for high-frequency modulations of the time delay 
for different values of the modulation amplitude. The 
planes of the control domains are parametrized by the 
feedback gain and the nominal delay of the control force 
for a fixed delay of the RDDE, and also, by the time de- 
lay of the original system and the nominal delay of the 
feedback control force for a fixed value of the feedback 
gain. The control domains are also determined for a low- 
frequency modulation in the plane of the feedback gain 
and the nominal delay of the control force for different 
values of the frequency of the modulation. The results 
show a significant enlargement of stability areas of VDFC 
with respect to TDFC within a certain range of the con- 
trol parameters, sometimes resulting in a complicated re- 
configuration depending on the type, the amplitude and 
the frequency of the delay modulation. The conclusions 
are summarized in Section IV. 



II. STABILITY ANALYSIS 

We consider a general nonlinear dynamical system de- 
scribed by a scalar autonomous RDDE in the form: 



x{t)^F[x{t),x{t-Ti)], 



(1) 



where Ti > is a constant delay time, and F is an arbi- 
trary nonlinear function of the state variable x, having a 
past dependence through the same state variable x but 
at Ti time units in the past. The presence of the delay 
term x{t — Ti) is a cause for the system ([T]) to be infinite 
dimensional, since a continuum of initial conditions over 
the time interval [—Ti, 0] is required in order to uniquely 



specify the future behavior of the system. The system 
possesses a set of fixed points {x*} that are solution of 
F[x*{t),x*{t — Ti)] — 0, and the stability of a particular 
fixed point x* can be obtained by linearizing Eq. ([T|) in 
the vicinity of x* . The linearized version of ([T]) around 
X* has a general form: 



x{t) ^ Ax{t) + Bx{t-Ti), 



(2) 



where A and B are real constants. We made a coordinate 
transformation from x to a; according to x{t) — x{t) — x* 
such that the fixed point is at the origin as expressed in 
the new coordinate. Employing the usual ansatz x{t) '^ 
exp(Ai) in ([2]) we obtain the characteristic equation: 



X^ A + Be 



-ATi 



(3) 



This is a transcendental equation in A, possessing a 
countable infinite set of complex solutions {Xi} defining 
the eigenvalues of the fixed point at the origin. The ori- 
gin is stable if and only if each A^ has a negative real 
part, it is unstable if at least one Xi has a positive real 
part, and it is marginally unstable if the largest real part 
of all the eigenvalues {Xi} is zero. 

The goal of this paper is to investigate the possibility of 
stabilization of the unstable fixed point x* of the system 
dD) by applying a Pyragas-type feedback force u{t) with 
a variable time delay [37| : 



u{t) 
r{t) 



K[x{t-T{t)) 
T2+ef{iyt), 



<t% 



(4) 
(5) 



such that for a given set of control parameters 
{K,T2,s,v} the unstable fixed point x* of the unper- 
turbed system ([ij becomes stable in the presence of the 
feedback term ^. The control parameter K is the feed- 
back gain characterizing the strength of the feedback, 
and T(i) is the variable time delay. We will consider a 
variation in a form of a deterministic modulation around 
a nominal delay value described by the control parame- 
ter T2. We take the delay function / : M — )• [— 1, 1] to 
be periodic with zero mean, with e and v being the pa- 
rameters determining the amplitude and the frequency 
of the modulation, respectively. The form of the control 
force (HI)-® implies that since T{t) > 0, the values of 
the amplitude e are restricted to the interval [0,r2]. In 
the presence of the control force ((4]), the system ([T]) has 
the form: 



x{t) =F[x{t),x{t^Ti)]+u{t), 



(6) 



and the linearized version around a;* in terms of the new 
coordinate x is: 



where 



x{t) = Ax{t)+Bx{t-Ti) + u{t), 



u{t) = K[xit-T{t)) -x{t)]. 



(7) 



(8) 



The stability of the origin can be inferred by numeri- 
cally integrating the linear variable-delay system ©-([H]) 
for different values of K^ T2, e and v, thus determining 
the domains in the {K, T2, s, v) hyperspace for which the 
stabilization becomes possible. 

For a sufficiently large variation of the time delay r(i), 
the stability of the linear variable-delay system ([T])-® 
becomes amenable for analytical treatment j46j |. From 
the stability point of view, if the frequency of the delay 
variation v is sufficiently large, then the linear system ([7]) 
with a variable time-delay (jHl) behaves as the following 
time-invariant system with a distributed delay (Theorem 
Al, Appendix A): 

x{t) = Ax{t) + Bx{t~Ti) + 



w{t]) x{er] + t-T2)dr]- x{t) , (9) 



K 



with w being the weight related to the probability dis- 
tribution of the delay function / in the interval of its 
periodicity, satisfying /_^ w{rii) drj ^ 1 (see Table I). The 
stability of the distributed delay system (jH]) is determined 
by the roots A^ of its characteristic equation: 



X^A + Be 



-ATi 



K 



'XT2 



5(Ae) - 1] , 



(10) 



C is a smooth complex function defined 
g{Xe) = [ w{r^)e^'^dv. (11) 



where g : 
as: 



In this sense, the solutions {Ai} determining the stabil- 
ity of the comparison system ([9]) can be considered as ef- 
fective eigenvalues describing the overall stability of the 
original variable delay system dZ])-®, providing that the 
delay frequency f is large compared to the system's dy- 
namics. Numerical simulations showed that the thresh- 
old for the frequency i^ above which this type of compar- 
ative analysis becomes valid needs not to be very high, 
and that its value depends on the actual system under 
investigation. 



A. Stability of the unperturbed system 

In the absence of control, the stability of the fixed point 
X* is determined by the roots of the characteristic equa- 
tion ([3]). Let Hq be a function of A defined as: 



Ho{X) = X-A- Be 



-ATi 



(12) 



With the aid of this characteristic quasipolinomial Ho{X), 
Eq. ([3]) can be written as Ho{X) = 0. We would like to 
find the range of the values for A, B and Ti for which x* 
is stable. 

Since Hq(X) is a smooth function on A, it is useful to 
consider the behavior of Hq(X) as A changes continuously 



over the real interval [0, - 
this interval, we have: 



-oo). Specifically, at the ends of 



lim Hq{X) = +00, 

A— >oo 

lim Ho(X) = -(A- 

A^0+ 



B). 



(13) 
(14) 



li A + B > 0, then Hq changes its sign at least once as A 
sweeps along the positive real axis. Consequently, there 
exists at least one positive real root of the characteristic 
equation Hq(X) = 0, rendering the fixed point unstable 
for any Ti. li A + B = 0, then A = is a root of the 
characteristic equation ®, and the fixed point is unsta- 
ble, or at least marginally unstable. Hence, a necessary 
(but not sufficient!) condition for stability of the fixed 
point is: 



A + B <0. 



(15) 



Taking into account that the boundary between stabil- 
ity and instability (the threshold of control) occurs when 
the maximal value from all the real parts in the set of 
solutions {Xi} is zero, we look for a solution of Eq. ([3]) in 
the form X = iuj, uj € M., and separate real and imaginary 
parts of the resulting equation to obtain: 



A = Bcos(wTi), 
-uj = Bsm{ojTi). 



(16) 

(17) 



[We stress that a zero on the imaginary axis for some set 
of parameters A, B and Ti does not necessarily mean that 
all the other zeros of the characteristic polinomial Ho{X) 
for the same set of parameters have negative real parts. 
The stability boundary is just one set of solutions of Eqs. 
(rT6|) - p7)) .] By eliminating the trigonometric terms from 
the last pair of equations, we get: 



B^ -A' 



(18) 



from which we conclude that Eq. ([3]) can have a solution 
for A on the imaginary axis if and only if |i3| > \A\. Tak- 
ing into account that Ti > 0, from Eq. (|T6)) we obtain: 



Ti 



Atccos{-A/B) + 2mr 

VS2 - ^2 ' 



(19) 



where n is a nonnegative integer, and Arccos denotes the 
principal value of the arccosine function. Obviously, the 
first value of Ti for which Ho{X) = has a solution for A 
on the imaginary axis is: 



Ti* = 



Arccos (- A/ S) 



(20) 



which follows from Eq. ([T^ by setting n = 0. The 
behavior of the real part of A at the values for Ti in Eq. 
P^ is determined by the derivative dX/dTi at X — ioj. 
By implicit differentiation of Eq. ([3]) with respect to Ti , 
we obtain: 



dX 
dJ\ 



XBt 



-ATi 



A(A - A) 



1 + BTie-^^i l + Ti{X-A)' 



(21) 



TABLE I: A representation of the delay function /, the weight w of the distributed delay system, and the function g, corre- 
sponding to three different types of delay modulations. By Iq we denote the modified Bessel function of the first kind of order 
zero, Jo is the Bessel function of the first kind of order zero, and S is the Dirac delta function. 



Type 


fit) 


w{t) 


ff(Ae) 


g{iuj£) 


Sawtooth wave 




1 
2 


sinh(Ae) 
Ae 


sin(a;e) 




sin(f) 


i 


/o(Ae) 


Jo (i^e) 




TtVI - t2 


Square wave 


f -1, te[o,Tv) 

\ 1, te [n, 2tt) 


S{t-l) + 5{t+l) 
2 


cosh(Ae) 


cos(aje) 



from which at X — iio we get: 



Re 



dX 
dTi 



U! 



A— icj 



(1 - ATi)2 + (wTi) 



(22) 



Since the sign of this derivative is always positive, the 
sign of the real part of A switches from negative to pos- 
itive when the zero of the characteristic quasipolinomial 
Ho{X) crosses the imaginary axis. On the other hand, 
as an implication of the Rouche theorem, the number of 
roots (counting multiplicity) on the complex right half 
plane (RHP) and the number of roots on the complex 
left half plane (LHP) can be changed (or, more correctly, 
interchanged) o nly if a zero appears on or crosses the 
imaginary axis |47l. l48l| . As a consequence, in the case 
under consideration \B\ > \A\, all the zeros of the char- 
acteristic quasipolinomial Ho{X) lie on the LHP if Ti is 
in the interval [0,T-^) providing that all the zeros were 
on the LHP before the first crossing of the imaginary 
axis has occured. However, this is evidently not true for 
other intervals separated by the corresponding values of 
Ti given by Eq. ((T9)) for n > 0, since the first zero- 
crossing of the imaginary axis occurs for Ti = Tj", and 
according to Eq. (|22|) every crossing is from the LHP to 
the RHP. 

In the case A> 0, the necessary condition for the sta- 
bility of the fixed point x* is B < —A [see Eq. (fTS])]. 
which is an interval of B that belongs to the range 
\B\ > \A\ for which the characteristic equation (jH]) can 
have a solution on the imaginary axis. From the pre- 
vious discussion, the possibility for all the zeros of the 
quasipolinomial Ho{X) to lie on the LHP necessary imply 
Ti e [0,Tj*). Since for Ti = the characteristic equation 
([3]) is reduced toA = A-|-i?<0, and since the crossing of 
the imaginary axis occurs for Ti = Tj* , we conclude that 
all the zeros {A^} have negative real parts in this case if 
and only ii B < ~A and Ti e [0, T*). 



In the case A < 0, the necessary condition for the 
stability of the fixed point is i? < jAj. In the subinterval 
B S [— 1^1, 1^1), the characteristic quasipolinomial p^ 
cannot have a zero on the imaginary axis. Choosing B = 
0, from ([3|) we have X = A < 0. Since crossing of the 
imaginary axis does not occur for this subinterval of B, 
it follows that aU the zeros {AJ for B e [-|A|, |A|) lie 
on the LHP for any Ti > 0. On the other hand, in the 
range B < —\A\ the characteristic quasipolinomial (fT2|) 
can have a zero on the imaginary axis. Putting Ti = 
in ([3]) we obtain X ^ A + B < 0, which means that when 
B < —\A\ all the zeros {Xi} have negative real parts when 
TiG[0,Ti*). 

The results are summarized with the following theo- 
rem: 



Theorem 1 Let the linear RDDE: 



x{t) ^ Ax{t) + Bx{t-Ti), 



A,B£ 



be a result of linearization of a corresponding nonlinear 
RDDE with a constant delay Ti : 

x{t) ^F[x{t),x{t-Ti)] 

around some fixed point x* of the latter expressed in co- 
ordinates in which the fixed point is at the origin. Fur- 
thermore, let Tj* > &e a real positive constant defined 
as: 



TI 



Arccos(-A/B) 

VS2-A2~' 



Then, the fixed point x* is locally asymptotically stable 
in each of the following cases: 

(a). B <-\A\ and Tie [0,r*); 



(b). B e [-\Al \A\), A<0 andTi> 0. 
Otherwise, x* is unstable. 



B. Stability under variable-delay feedback control 
(high-frequency modulation) 

In the following, we consider the modulation frequency 
v to be above the threshold, allowing an analysis of the 
variable delay system dZ])^® as a distributed delay sys- 
tem (|n]). When the control is switched on, the stability 
of the fixed point x* is determined by the roots {A^} of 
the characteristic equation pUj) . If we define: 

H,{X) = X-A~Be-^'^' +K [l~ e-^'^'-g{Xe)] , (23) 

then Eq. PH)) can be rewritten as H^{X) = 0. Assum- 
ing that in the absence of control, the parameters A, B 
and Ti of the unperturbed system are such that x* is un- 
stable, we look for the values of the control parameters 
K , T2 and e for which the fixed point is stabilized. In 
other words, we would like to find the set of points (i. 
e. to determine the domain of control) in the parameter 
space {K, T2, e) for which all the zeros of the characteris- 
tic quasipolynomial -ffe(A) lie on the LHP, while, at the 
same time, the characteristic quasipolynomial Ho{X) of 
the unperturbed system has at least one zero in the RHP. 
Before we proceed with the analytical description of 
the control boundaries, it is interesting to consider the 
behavior of H^{X) as A changes continuously over the 
positive real axis. Taking into account that g{0) — 
/_! w{v) drj — 1, from Eq. (j23p we obtain: 



lim -ffe(A) = +00, 

A— voo 

lim HJX) = -(A + B), 

A-J-0+ 



(24) 
(25) 



which coincide with the limits (IT^ - p^ for the char- 
acteristic polynomial Ho(X) of the unperturbed system, 
leading to the same necessary condition ([15]) for stabil- 
ity of the fixed point. Since (fT5|) does not include the 
dependence on the control parameters K, T2 and e, we 
conclude that VDFC is unsuccessful for any values of the 
control parameters if the linearized version ([2]) of the un- 
perturbed system around x* is such that A-\-B > 0. This 
important result is expressed in the following theorem. 

Theorem 2 Let x(t) = Ax{t) + Bx{t-Ti), A, B e M, 
be a linearization around the fixed point x* of the corre- 
sponding nonlinear RDDE with a constant delay Ti . If 
A-\- B > Q, then the variable- delay feedback control ^- 
(0) cannot stabilize the unstable fixed point x* for any 
value of the control parameters K , T2 and e. 

The limitation of the VDFC method imposed by Theo- 
rem 2 is a kind of an analogue to the odd-number lim- 
itation |30l-l32| in the case of delayed feedback control 
of systems described by ordinary differential equations. 



whose validity was recently refuted [33|-|36j for the case 
of unstable periodic orbits. 

Exact analytical description of the domains of success- 
ful control in the parameter space (K,T2,e) is difficult 
for the characteristic Eq. pO)) due to the complexity 
of the terms involving the dependence on A. Thus, one 
should solve Eq. ([T0| numerically in order to calculate 
the control domains. To this extend, it is possible to ob- 
tain expressions for the parametric representation of the 
control boundaries parametrized by a Hopf frequency uj. 
Substituting A = iw in Eq. ([TU)) and separating real and 
imaginary parts, we obtain: 

Kg{iuj£) COSUJT2 = K - A- BcosujTi, (26) 
Kg{iuje) sin LUT2 = —u — BsinujTi. (27) 

Elimination of T2 from the last pair of equation yields a 
quadratic equation in K: 

[1 - giiujef] K^ -2{A + B cos wTi) K 

+ {A + BcosujTif + {uj + B sin wTif =0, (28) 

which can be solved for K in terms of uj to get: 



K{uj) 



A + BcosujTi 



± 



1 



1 ~ [g{iuje)] 1 - [g{iuje)] 

[g{iuje)f {A + B cosujTif 



+ {[g{iuje)f -l){uj + BsinuTiy' 



1/2 



(29) 



On the other hand, by dividing (|27|) and ([26| . we obtain 



T2{u;) = - 

UJ 



Arctan 



—UJ — B sinujTi 



K- A 



± niT 



BcosujTi 

" (30) 

which, together with Eq. ((29| . describe the stability 
boundary for a fixed e in the {K, T2) plane, parametrized 
by w. 

It is also useful to study the stability boundaries of 
the controlled system for a fixed feedback strength K 
and modulation amplitude e in the parameter plane of 
the two delay times (Ti,T2). Following the idea in Ref 
[ill, we rewrite Eq. dTU]) as: 



l + a(A)e~^^i +6(A)e 
where a(A) and 6(A) are given by: 

B 



-XT2 



0, 



a(A) = 
6(A) = 



A-K-X' 

KgjXe) 
A-K-X' 



(31) 

(32) 
(33) 



At the control boundary (A = iuj) the three terms in Eq. 
([31]) can be considered as three vectors in the complex 
plane, with the corresponding magnitudes 1, 10(1^)! and 
|6(ia;)|. According to Eq. (I5T|) . the sum of these vectors 
is a zero vector, thus forming the triangle shown in Fig. 




Ser 



FIG. 1: Diagram in the complex plane 2:, associated with the 
derivation of the parametric representation of the stability 
boundary in (Ti,T2) plane. 



1. From Fig. 1, it is straightforward to obtain the para- 
metric representation of Ti and T2 on the Hopf frequency 

y^(^) ^ Arg [a{iu)] + {2u ^ l)7r ± ^1 ^ ^ 



T2(w) 



u — Uq,Uq + 1,mJ + 2 . 



Arg [b{iuj] + {2v - 1)tt t 6>2 
v = Vq,Vq +1,Vq +2..., 



>0, 



(34) 



(35) 



terms of lo: 










j^i ^ A + Bgiiue) , 1 

l-[g(ia.e)f 1 - [5(*^e)f 


{A + Bg{iuje)f 


+ i[g{iue)]^ - l)iA^ + B^ - u^) 


1/2 

(41) 


T{u) = - 

UJ 


\. f -^ \ 1 


(42) 


r^^"^^\K-Aj '""J 



When e = 0, VDFC reduces to the usual Pyragas con- 
trol scheme (TDFC) with a constant delay T2. Since 
TDFC is a special case of VDFC when the modulation 
of the control delay in the feedback force is absent, the 
parametric representations of the control boundaries for 
TDFC simply follow from the ones derived in the case of 
VDFC by letting e = (or, equivalently, g{0) = 1) in the 
corresponding equations. For example, from Eqs. (j26p - 
(j27p with e = we obtain the parametric representation 
of the TDFC boundary in the {K, T2) plane parametrized 
by w: 



K{u) 
T2{lo) = - 



[A + BcosuTif + (lo + BsmojTi) 



Arctan 



2{A + Bcosu}Ti) 

—Lu — B sin ujTi 

K -A- BcosujTi 



± nn 



,(43) 
.(44) 



It is interesting to note that when Ti — T2 ^ T in the 
case of TDFC, the corresponding characteristic equation 
can be written as: 



where Uq and Vq are the smallest possible integers such 
that the corresponding values of Ti and T2 are all non- 
negative, and 6*1,^2 G [0,^] are the internal angles of 
the triangle shown in Fig. 1 calculated from the law of 
cosines as: 



X = A' + B', 



-XT 



di = Arccos 
02 = Arccos 



l + |a(zcj)p-|b(zaj)P' 

2\a{iuj)\ 
l + |b(icj)p-|a(ic^)P 

2\b{iuj)\ 



(36) 

(37) 



In the case when the nominal delay T2 of the feedback 
control force coincides with the delay of the original sys- 
tem Ti, the characteristic Eq. ([T0|) is reduced to: 



\- A + K~[B + Kg{\e)] e 



-XT 



0, 



(38) 



where we use T — Ti — T2- At the stability boundary 
(A = ibj) the last complex equation can be represented 
as a pair of two real equations: 

[B->rKg{iuje)]cos.ujT = K - A, (39) 

[B + Kg{iuje)]smujT = -uj, (40) 

which can be manipulated to obtain a parametric repre- 
sentation of the control boundary in the {K, T) plane in 



(45) 



where A' = A — K and B' = B + K. Noting the equiva- 
lency between Eq. (j45|) and Eq. ((3|) , the exact analytical 
description of the stability domain in this case immedi- 
ately follows from Theorem 1. 



III. NUMERICAL EXAMPLE 

To test the VDFC method for stabilization of unsta- 
ble steady states in chaotic RDDE systems, we will use 
the paradigmatic Mackey-Glass system introduced as a 
model for regeneration of blood cells in patients with 
leukemia |49l - l53J |. The Mackey-Glass equation in the 
presence of VDFC states: 



x{t) 



ax{t — Ti) 
l + [x{t-Ti)]' 



bx{t)+u{t), (46) 



where u{t) is given by Eqs. ©-([S]). Here x{t) is a 
concentration of circulating blood cells, and a, 6 and 
c are parameters of the free running system, involved 
in the description of the dependence of the produc- 
tion/destruction of the blood cells as a function of x{t) 
and x{t — Ti), respectively. We will consider the typical 
values a = 0.2, h = 0.1 and c = 10. 
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FIG. 2: Representative samples of phase plots x{t) vs. 
x{t — Ti) for the uncontrolled Mackey- Glass system at dif- 
ferent values of Ti : (a) Ti = 4 - the trajectory is attracted to 
the stable equilibrium point x'^ = —1; (b) Ti = 8 - the trajec- 
tory approaches a limit cycle; (c) Ti = 15 - the attractor has 
evolved into a period-2 cycle; (d) Ti = 23 - chaos. The simu- 
lations were performed using the MATLAB routine dde23 for 
integrating delay-differential equations with constant delays. 



In the absence of control [u{t)—0], the system (|46)) has 
a set of three fixed points xl = 0, Xj = +1 and Xj = —1 
being solutions of: 



ax 



1 + x* 



hx* =0. 



(47) 



The stability of each x* is obtained by linearizing the 
unperturbed system around x* , leading to Eq. ([2]) with: 



A 



B ^a 



l-h(l-c)a;* 



(48) 



and the corresponding characteristic equation is given by 
Eq. ©. For x*^ = 0, we have A = -b = -0.1 and B = 
a — 0.2. Using Theorem 1 we deduce that the fixed point 



is unstable for any Ti. For x 



2,3 



±1, we have A 



-b = -0.1 and B = a(2 - c)/4 = -0.4, indicating that 
this pair of fixed points are characterized by the same 
type of stability. From Theorem 1 we conclude that X2 3 
are stable if and only if Ti e [0,4.7082). Figure 2 shows 
the trajectory of the unperturbed system in x(t) vs. x{t— 
Ti) coordinate space for four different values of Ti. Panel 
(a) shows the evolution of the system for Ti — 4. Since 
for this value of Ti the fixed points Xj 3 are stable, the 
preference of the system towards 0:2 — +1 ^^ ^3 — ~1 
depends on the initial conditions. Panels (b)-(d) in Fig. 
2 correspond to Ti = 8, 15 and 23, respectively, showing 
the growth of the limit cycle through a period-doubling 
bifurcation sequence, and the eventual appearance of a 
chaotic attractor. 



In performing the stability analysis under VDFC, we 
will first consider a high-frequency modulation of the con- 
trol delay r(t). The limitation imposed by Theorem 2 
asserts that the fixed point xl cannot be stabilized with 
VDFC for any values of the control parameters K, T2 
and e. The validity of this assertion has been verified 
by the numerical simulations, showing the absence of the 
domains of successful control in the corresponding para- 
metric planes. On the other hand, the stability of the 
fixed points X2 ^ — ±1 is determined by the roots {Xi} of 
the characteristic Eq. ([TU| with A and B given by Eq. 
(|48p . Even though there exists an infinite number of roots 
Xi of Eq. (|10p . only a finite number of them have real 
parts greater than a given constant. A computation of 
the rightmost characteristic roots with large enough ac- 
curacy is a nontrivial nonlinear eigenvalue problem, and 
there exist several effective methods to compute this part 
of the spectrum, e.g. by a discretization of either the 
time integration operator or the infinitesimal generator 
associated with the delay system |54| - [58{ . Since the sta- 
bility properties of the controlled system are determined 
by the characteristic roots with the leading real part, it 
is enough to employ a simple root-finding numerical al- 
gorithm based on the Newton-Raphson iteration method 
with a suitable chosen grid of starting values. For this 
purpose, we first make an implicit plot of the real and 
the imaginary parts of the characteristic Eq. pH)) in the 
complex A plane to visualize the approximate location of 
the roots as intersecting points between the correspond- 
ing curves. In this way we obtain a coarse estimate of 
the location of the rightmost eigenvalues, the knowledge 
of which is then used to choose an appropriate grid of 
starting values encompassing this location. 



By numerically solving Eq. (1101) with the aforemen- 
tioned procedure, we obtain the domains of successful 
control in the parameter plane (iiT, T2) for a fixed delay 
Ti and for different values of the amplitude e. The results 
are shown in Figs. 3 and 4. In the numerical calculations, 
we choose Ti — 23 for which the original system is in a 
chaotic regime (see panel (d) in Fig. 2), having a positive 
value of the largest Lyapunov exponent (LLE = 0.00973) 
[59| . The shaded areas (color online) correspond to 
the set of control parameters {K,T2) for which the max- 
imum of the real part of the characteristic eigenvalues 
{Xi} is negative (max[Re{Ai}] < 0), rendering the control 
successful. The values of max[Re{Ai}] are given by the 
grayscale (colorscale online) on the right in each figure, 
and the control is more robust as max[Re{Ai}] is more 
negative. The stability islands are surrounded by a "sea" 
of instability represented by the white region, for which 
the real part of the leading characteristic eigenvalue is 
positive (max[Re{Ai}] > 0). The "coastline" between 
stability and instability (the stability border) is given in 
a parametric form via Eqs. (HHl-dSni) for e > (VDFC), 
and via Eqs. gSj-dMl) for e = (TDFC). Panels (a) 
through (d) of Fig. 3 correspond to the modulation of 
the feedback delay T{t) in a form of a sawtooth-wave, 
with amplitude values e = 0, 0.5, 1 and 2, respectively. 
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FIG. 3: (Color online) Domains of successful VDFC control 
in the [K,T2) plane for the unstable equilibria ^2,3 = ±1 
in the chaotic Mackey-Glass system (Ti = 23). The control 
delay T(t) is modulated with a sawtooth- wave, and the values 
of the modulation amplitudes are: (a) £ = (TDFC); (b) 
e — 0.5, (c) e — 1; (d) e = 2. Combinations of K and T2 
where VDFC successfully stabilizes the fixed points 3:2,3 = ±1 
are plotted in graytones (colortones online). Note the shifts 
of the origin along the T2 axes by an amount equal to e due 
to the limitation T2 > e. 



Panel (a) reveals the structure of the stability domain for 
e = (TDFC). For the current choice of Ti, and also in 
general, there exists a stability region for relatively small 
Ti with a complex structure, and a resonance island en- 
compassing T2 = Ti = 23 for which the control is most 
robust and can be achieved with smaller values of K. As 
e becomes larger than zero (VDFC, panels (b)-(d)), the 
structure of the stability domain is reconfigured, result- 
ing in a significant enlargement of the area of successful 
control. This enlargement is also observed for other de- 
lay modulations. In panels (a)-(b) of Fig. 4 we show 
the calculated stability domains for a sine-wave modu- 
lation for e = 1 and 2, respectively, and (c)-(d) are the 
corresponding panels for a square-wave modulation. We 
note that for larger values of e in the case of a square- 
wave modulation, the stability area eventually spreads 
into several clearly distinguished stability islands, whose 
position is changing in an oscillatory manner as e further 
increases. 

In Fig. 5 we show the stability domains in (Ti,T2) 
plane, fixing the feedback gain value a.t K = 0.5. Panel 
(a) depicts the case when the modulation is absent 
(TDFC, e = 0), and panels (b)-(d) are related to saw- 
tooth, sine and square-wave modulations, respectively, 
with e — 2. The diagrams show the typical enlargement 
of the stability area for VDFC with respect to TDFC. 
The parametric representation of the stability boundary 
is given by Eqs. (|54 l) - (f35l) . 

To verify the analysis in the previous paragraphs, we 






I 0.00 
-0.05 

I -0.10 

-0.15 

- -0.20 



10 20 30 40 



10 20 30 40 





10 20 30 40 

T2 T2 

FIG. 4: (Color online) (a), (b) Stability domains of a;2_3 — ±1 
in the {K,T2) plane for the VDFC-controUed Mackey-Glass 
system with Ti = 23. The delay modulation is in a form 
of a sine- wave with e = 1 (panel a) and e — 2 (panel b). 
(c), (d) Corresponding stability domains for a square- wave 
modulation. Note the shifts of the origin along the T2 axes 
by an amount equal to e due to the limitation T2 > e. 
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FIG. 5; (Color online) Domains of successful control in the 
{Ti,T2) plane for the unstable fixed points 12,3 = il in the 
Mackey-Glass system. The feedback gain is fixed at K = 0.5. 
(a) Stability diagram for e = (TDFC). (b)-(d) Respective 
stability diagrams for sawtooth, sine and square- wave modu- 
lations with e — 2 (VDFC). Note that the minimum value of 
the T2-axis in panels (b)-(d) is T2 = 2 due to the limitation 
T2 >£. 



performed a computer simulation of VDFC for the fixed 
points 2:2 3 — ±1 by numerically integrating the sys- 
tem p6)) for different delay modulations. The results are 
shown in Fig. 6. Panels (a), (c) and (e) depict the dy- 
namics of the variable x{t) for sawtooth, sine and square- 
wave modulations, respectively, and panels (b), (d) and 
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FIG. 6: VDFC applied to the chaotic Mackey-Glass system using difTerent modulations of the delay-time r(f). The parameters 
of the uncontrolled system are: a — 0.2, b — 0.1, c = 10, Ti = 23. (a), (b) Time plots of the variable x{t) and the feedback 
signal u{t) for a sawtooth-wave modulation, indicating a successful control of the unstable fixed point at X2 = -1-1. (c), (d) 
Stabilization of the unstable equilibrium ai x^ — —1 with a sine- wave modulation, (e), (f) Time-series for a square- wave 
modulation stabilizing the unstable point at Xg — —1. In each case, the control parameters were: K = 2, e = 2, T2 = 18 and 
1^ — 5. The control was activated at t = 500. The total time span shown in each panel is 1500 time units. The simulations 
were performed using the MATLAB routine ddesd for integrating delay-differential equations with general delays. 



(f) show the corresponding time-series of the feedback 
signal u{t). In each case, the control parameters were 
chosen as -fC = 2, e = 2, T2 = 18 and v = 5, fixing the 
delay of the uncontrolled system at Ti = 23 for which 
the system is chaotic. We note that for these parameter 
values, the control via TDFC (e = 0) is unsuccessful for 
any K, as can be perceived from the stability domain de- 
picted in panel (a) of Fig. 3. Also, since Xj 3 have identi- 
cal set of characteristic eigenvalues, they share common 
domains of successful control. However, they have dif- 
ferent basins of attraction, and the preference of control 
towards either ccj or Xg depends on the initial conditions. 
In panels (b) , (d) and (f ) we see that the feedback signal 
u{t) vanishes when the stabilization of the fixed point is 
achieved, suggesting noninvasiveness of VDFC, which is 
a consequence of the form of the control force in Eq. ^ , 
since x(t — T{t)) = x(t) if the fixed point is stabilized. 

When the frequency z/ of the delay modulation is below 
the threshold (low- frequency modulation), the approxi- 
mation of the variable-delay system with a distributed- 
delay system is not covered by Theorem Al, and, hence, 
the control domains cannot be calculated from the char- 
acteristic Eq. ([To]). However, the stability domains in 
this case can be obtained by numerically integrating the 
linear variable-delay system ([T])-® for different values 
of the corresponding control parameters. In Fig. 7 we 



show the results of such a simulation in the paramet- 
ric plane {K,T2) for Ti — 23 and e — 2, taking the 
time-modulation of T{t) in a form of a sawtooth wave. 
Different panels of the figure correspond to different val- 
ues of the delay frequency v. (a) i^ = 1.4, (b) v = 1.8, 
(c) 1^ = 1.9, (d) 1^ = 2.0, (e) 1^ = 2.2, (f) j. = 2.4, (g) 
1/ — 2.6, (h) ly — 3.0. The combinations (if, T2) leading 
to a successful stabilization of the unstable fixed points 
X2 3 = ±1 are marked in black. It is observed that when 
the modulation frequency is about ly = 3.0 (panel (h)), 
the structure of the stability domain fairly resembles the 
stability domain for a high-i^ modulations obtained from 
the characteristic Eq. PH)) (compare with panel (d) in 
Fig. 3, noting the different scales on the T2 axis). As 
expected, the simulations show that this resemblance be- 
comes improved as v attains higher values. On the other 
hand, the structure of the stability domain is gradually 
changing as 1^ becomes smaller than 1^ — 3.0 (panels (a)- 
(g)), resulting in a reconstruction of the main domain and 
a birth of many small stability islands, clearly notable 
for larger nominal delays T2 and approximately centered 
about those T2 which are odd multiples of n/v. The 
emergence of this additional domain structure could be 
due to a resonance between the delay frequency v and the 
intrinsic frequencies of the uncontrolled system, which 
are infinite in number. The distance between these res- 
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onance islands (« 2tt/v) becomes wider as v decreases, 
and they become less pronounced for lower values of v. 
It can be noticed that the appearance of these resonance 
islands allows stabilization of the unstable equilibria for 
much larger nominal delays T2 in comparison to the val- 
ues of T2 for a high-u modulation. The simulations show 
that the range of the delay frequency parameter con- 
taining these resonance islands is strongly dependent on 
the system parameters (e.g., the modulation amplitude) 
and on the type of the delay modulation, and that this 
range of v may not be continuous as in the current case, 
but it may consist of several different subintervals spread 
throughout the entire i/-interval below some sufficiently 
high frequency and encompassing some of the values of v 
coinciding with the eigenfrequencies of the uncontrolled 
system (see Fig. 8). 

To check if the limitation of the control method as- 
serted by Theorem 2 remains valid for low-frequency 
modulations, we have performed numerical simulations 
to determine the domains of successful VDFC control of 
the unstable equilibrium x\ — 0, which has been shown 
uncontrollable via Theorem 2 for high-frequency modu- 
lations. The simulations in this case show the absence 
of the control domains in the corresponding parametric 
planes, suggesting the validity of Theorem 2 in the entire 
frequency range. 



Moreover, numerical simulations suggest that this limi- 
tation is also valid for low-frequency modulations. Nev- 
ertheless, in lack of any analytical tool to treat a low- 
frequency modulated VDFC, any general statement con- 
cerning the generalization of Theorem 2 to the whole fre- 
quency range should be taken cautiously, as well as the 
related observations concerning the aforementioned res- 
onance phenomenon. 

Putting the observations related to low-frequency mod- 
ulation of the control delay on a firm mathematical ba- 
sis constitutes an interesting subject for a future study. 
Other possible directions for future consideration would 
be stabilization of unstable steady states by VDFC in 
other types of DDE systems (e.g. systems described by 
neutral delay-differential equations 43*), in systems de- 
scribed by partial differential equations, and, also, imple- 
mentation of the control method to stabilize unstable pe- 
riodic orbits by a suitable choice of the delay modulation 
in order for the control method to stay noninvasive. An 
example for such a modulation in the latter case would 
be a periodic change of the control delay between T and 
2T where T is the period of the orbit to be stabilized 
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IV. CONCLUSION 

In conclusion, we have shown that variable-delay feed- 
back control allows stabilization of unstable steady states 
in a class of scalar retarded time-delayed systems, repre- 
sented by the chaotic Mackey- Glass system, over much 
larger domain of parameters in comparison to the usual 
Pyragas' delayed feedback control scheme. The analy- 
sis showed that the enlargement of the control domain 
may undergo a complex rearangement depending on the 
type and the frequency of the delay modulation. It is no- 
ticed that the enlargement of the control domain for high- 
frequency modulation of the delay is more pronounced 
when the variable delay is a continuous function of time 
in contrast to the case of variable delay function with a 
discontinuity leading to complex stability domain struc- 
ture of a lesser magnitude. In the case of low-frequency 
modulation of the delay, we notice a complex rearrange- 
ment of the control domain, resulting in an appearance of 
extra stability islands, probably a consequence of a res- 
onance between the frequency of the variable delay and 
the eigenfrequencies of the uncontrolled system. This 
resonance effect allows successful stabilization of the un- 
stable fixed point for much larger nominal delays with 
respect to the situation when the frequency of the delay 
variation is above the threshold. 

Limitation imposed by Theorem 2 shows that VDFC 
method fails to control certain unstable steady states 
for any value of the feedback control parameters in the 
case when the frequency of the delay modulation is high. 
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the discussion related to the resonance phenomenon in 
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discussions with B. Fiedler on the limitations of the "odd- 
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Appendix A 

The stability of a linear (or linearized) RDDE system 
with a fast-varying delay can be obtained by studying 
the roots of the characteristic equation of the related 
time-invariant distributed delay system. The correctness 
of this approach is guaranteed only if the frequency of 
variation of the delay is large compared to the system's 
dynamics. A precise formulation of these assertions con- 
stitutes the following theorem: 

Theorem A.l Consider the linear system of variable 
delay differential equations: 



— x(i) = A-x(t)+B-x(<-ri) + u(i), (Al) 
at 



u(i) - K.(x(i-T(t))-x(i)), 
T{t) = T2 + ef{vt), 



(A2) 
(A3) 



where A,B,K e R^^^ are constant matrices, x(t) e 
M"'^^"'^, and f : M ^- [—1,1] is a periodic function with 
zero mean and period 2n, max/ — 1, and min/ = — 1. 
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FIG. 7: The stability domains in {K,T2) plane related to the unstable steady states X2^z ~ il of ^^e chaotic Mackey- Glass 
system (T\ — 23) for a low frequency modulation of the time delay r(t). The delay modulation is in a form of a sawtooth- wave 
with e — 2, and the value of the modulation frequency is: (a) f = 1.4, (b) v = 1.8, (c) v = 1.9, (d) v = 2.0, (e) u = 2.2, 
(f) v — 2.4, (g) v — 2.6, (h) v = 3.0. The combinations [K,T2) leading to a successful stabilization of the unstable equilibria 
X2,3 = ±1 are marked in black. The characteristic eigenfrequencies of the uncontrolled system that lie in this interval of u 
are: 1.43, 1.71, 1.98, 2.25, 2.52 and 2.80. Note the appearance of the resonance islands at the right of the main structure as v 
becomes smaller than v — 3.0. Also note the shifts of the origin along the T2 axes by an amount equal to e due to the limitation 
T2 >£. 
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FIG. 8: A sample of the low-frequency control domains in {K, T2) plane corresponding to the unstable steady states 3:2,3 = il 
of the chaotic Mackey- Glass system (Ti = 23) for a square-wave modulation with e = 2. The value of the modulation frequency 
is: (a) V = 1.7, (b) v = 2.8, (c) v = 3.0, (d) v = 4.0, (e) v = 8.0, (f) v = 10.0, (g) v = 11.0, (h) v = 12.0. The range of v in 
which the resonance islands exist consists of several distinguished intervals encompassing some of the eigenfrequencies of the 
uncontrolled system: 1.71, 3.07, 7.99, 10.99. Note that although the modulation frequency v in panel (b) coincides with one of 
the eigenfrequencies of the uncontrolled system (~ 2.8), the resonance islands are hardly noticeable in this case. 
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Let e,T2,i' E Mq , and e < T2. Let the integrable function 
w : [— 1, 1] — > M'^ be defined by: 



Y{t)w{t)dt = 



27r 



a{f{t))dt (A4) 



for every continuous function a : [—1,1] — > K. // the 
comparison system: 



— xm = A • x(<) + B • x(t - Ti)- 
dt 



+ K 



r-t-Ts+e 



It-T-i-e 



j{{e-t + T2)/e)_ 



:{9) dO - x(i) 



(A5) 



is asymptotically stable, then the original system hAl\) - 
is globally uniformly asymptotically stable for large 



values of the frequency v of the modulation. 

Theorem Al is a restatement of the main result in Ref. 
f46| to accomodate the present discussion [60|, and its 
proof is based on an extension of the recently introduced 
trajectory-based proof technique 61]. According to The- 
orem Al, the stability of (jAll) under the variable-delay 
control force (|A2|) can be inferred from the stability of the 
analogous time-invariant system (jA5l) with a distributed 
delay, for sufficiently large values of the parameter v de- 
termining the frequency of the modulation. It is worth 
noting that Theorem Al can be generalized to include 
the most general case of multiple delayed feedback terms 
in the control force (jA2|) with different types of delay 
modulations [60] . The proof of this extension is straight- 
forward, following the lines of the proof given in [4a] . 

The comparison system (|A5|) can be recast in the form: 



dt 



x(i) = A • x(i) + B • x(i - Ti) + K 



w{r])x{er] + t — T2) drj — x(i) 



(A6) 



by making a change of the integration variable to the 
new variable rj through the relation 6 = erj + t — T2. 
Furthermore, by taking a{t) — 1 and a{t) = t in Eq. 
(|A4p respectively, we obtain the relations involving the 
weight function w: 



I w{t)dt = 1, 



tw{t)dt = 0. 



(A7) 
(A8) 



From Eq. (jA4|) . the weight w can be interpreted as the 
probability distribution of /(^), where ^ is uniformly dis- 
tributed over the interval [0,27r] (see Table I). 
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